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The Barrett-Crane spin foam model for quantum gravity provides an excellent setting for 
testing analytical and numerical tools used to probe spinfoam models. Here, we complete 
the report on the numerical analysis of the single 4-simplex area correlations begun in Phys. 
Lett. B670 (2009) 403-406, discussing the next-to-leading order corrections ("one-loop" cor- 
rections) with particular attention to their measure dependence, and the difference between 
the Gaussian and Bessel ansatzc for the boundary state. 



I. INTRODUCTION 

The current challenge for spinfoam models of quantum gravity [1] is the understanding of their 
semiclassical regime: how to recover general relativity and the subsequent quantum corrections 
due to the quantum fluctuations of the geometry (around the flat space-time) from the spinfoam 
amplitudes. A key proposal by Rovelli [2] for computing the graviton propagator through 
correlations of geometrical observables has provided a well-suited tool to probe the semi-classical 
behaviour of spinfoam models [3-8] (see also [9] for a review). 

The most studied spinfoam model until recently was the Barrett-Crane model [10, 11]. The 
corresponding spinfoam amplitudes are built from the lOj-symbol first defined in [10]. Its 
asymptotic limit, relevant to the semi-classical behaviour of the model, has been well-studied 
analytically [12-17]. The leading order is dominated by degenerate geometries [15-17] and the 
contribution corresponding to non-degenerate configurations is an oscillating sub-leading term of 
the asymptotic expansion. The key idea of Rovelli's proposal was to kill the degenerate leading 
order through destructive interference with the boundary state and then compute the "graviton 
propagator" using the proper non-degenerate oscillating term [2]. This was rigourously imple- 
mented in [4] where the graviton propagator was computed analytically for a single 4-simplex at 
leading order. These results were checked numerically in [7], which provided the first numerical 
confirmation of the sub-leading oscillatory behaviour of the lOj-symbol. There was a perfect 
fit with the analytical predictions of [4] at leading order in the asymptotic regime. Moreover 
the numerical simulations allowed us to probe the small scale behaviour of the Barrett-Crane 
model, showing that the correlations were regularized at the Planck scale. Indeed, although the 
correlations were found to decrease as expected as the inverse squared distance at large scale 
(above a hundred Planck lengths) , they reached a maximum at a few Planck lengths when going 
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to the zero-distance limit and then decreased again. This correlation peak provides a dynamical 
definition of a minimal distance and shows that the gravitational correlations are regularized at 
short scales (at least when looking at a single 4-simplex). 

Although the Barrett-Crane model gives this remarkable result, it does not lead to the right 
semi-classical behaviour. Indeed, although the correlations properly fall as 1/r 2 at large dis- 
tances, they do not have the correct tensorial structure [3, 5]. This opened the door to the 
new EPR-FK-LS class of spinfoam models which are better-behaved semi-classically and are 
based on a relaxed imposition of the so-called simplicity constraints of the Barrett-Crane model. 
These have been constructed algebraically [18, 19], with coherent states techniques [20-22], and 
using discrete path integral techniques [23, 24]. This new EPR-FK-LS spinfoam vertex (for 
the quantum 4-simplex) was shown analytically and numerically to propagate coherent wave 
packets of 3d geometry [25]. Its asymptotic behaviour improves on that of the Barrett-Crane 
vertex [26], and it is related to Regge calculus in area-angle variables (see [27]). Finally, the 
graviton propagator was computed in this framework for a single 4-simplex and was found to 
yield the correct tensorial structure at leading order [28]. Therefore, these EPR-FK-LS spinfoam 
models have a better semi-classical behaviour than the previous Barrett-Crane model and there 
is much hope that we will recover classical general relativity plus quantum gravity corrections 
in a continuum limit. 

On the other hand, the Barrett-Crane amplitudes still provide a sensible semiclassical dy- 
namics when restricted to area correlations on a single 4-simplex. And being much simpler than 
the EPR-FK-LS amplitudes, they provide a useful tool for numerical and analytical studies. 
Their simplicity comes from the absence of intertwiner degrees of freedom, which is ultimately 
also the reason why we can not look at more generic geometric correlations nor properly glue 
more 4-simplices together. The model is also the simplest non-top ological spinfoam model for 
4d quantum gravity, thus an arena for studying issues like semiclassical symmetries and renor- 
malization (under coarse-graining) in the formalism. 

In the present paper, we push further the numerical work started in [7]. In the next section, 
we start by reviewing the definition of the Barrett-Crane model and the ambiguity in defining 
the path integral measure. We also present a plot showing the degenerate contribution and 
two other sub-leading contributions which also mask the expected oscillations. Then, in Section 
III, we define the correlations. We restrict ourselves to factorized (homogeneous) 

boundary states. We recall the two types of ansatz for such states: the (phased) Gaussian 
ansatz originally introduced in [2] and the Bessel ansatz later defined in [4, 29]. The Bessel state 
is simply a refinement of the Gaussian ansatz allowing more precise analytical computations 
but the two states are asymptotically the same. We then present our numerical results, plots 
of the correlations in term of the area scale jo for different choices of the Gaussian width a and 
the parameter k labeling the ambiguity in the path integral measure (more precisely the face 
amplitude). We recover the typical behaviour already identified in [7]: a plot with a peak for 
some given small value of jo an d then the correlations going down as 1/jo for large jo- The 
precise leading order coefficient in front of 1/jo depends on a but is not sensitive to k while the 
precise position of the peak depends on both a and k. We also find a new oscillatory regime for 
intermediate values of jo, which is rather unexpected. 

In Section IV, we focus on the next-to-leading order (NLO) correction in the asymptotic 
regime, making use of new numerical methods [35] that allow us to obtain the required precision. 
We find a sub-leading term that goes like 1/jo- We estimate the coefficient for different choices 
of a and k and we conjecture a quadratic dependence on the measure ambiguity k. This NLO 
correction should correspond to a one-loop correction to the graviton propagator and Newton 
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potential (which defines the /3-function for Newton's constant at first order). In Section V, 
we discuss the small scale behaviour, interpreting the peak as defining a dynamical notion 
of a minimal distance. In Section VI, we turn to the Bessel state which introduces further 
interferences and enhances the oscillatory regime for small/intermediate values of j$. Finally, in 
Section VII, we discuss the relevance of our choice of boundary state and the limits of the power 
of numerical simulations for the study of spinfoam correlations. We end with an appendix which 
gives a summary of the new numerical techniques which make the NLO calculations possible. 



The object that we study is the area correlation on a single Riemannian 4-simplex, defined 
in [2] and further studied in [3, 4, 8]. This is a 10 by 10 matrix with components 



In the following we will assume that the reader has a certain familiarity with this quantity. For 
more details and motivation, as well as a description of its connection to the graviton propagator, 
we refer to the literature cited above. In a few words, the notation used in (1) is as follows: 
the quantities j\ are spins, i.e. non-negative half-integers, which label SU(2) representations. 
They are indexed by subscripts a, b and I which run over the 10 triangles of a 4-simplex. 
The normalization M is the same as the displayed sum, without the field insertions h(j a ) = 
(Aq — Aq) /Aq, which represent the fluctuation of the area A a = £p (2j a + 1) of the triangle a 
around its background value Aq = £p (2 jo + 1). The function ty q is a boundary state, which 
should be peaked on the geometry of a regular 4-simplex. We describe our choices of this state in 
Section III A. The model-dependent kernel K[j] encodes the dynamics. In this paper, we focus 
on the Barrett-Crane model [10], and summarize the properties of the kernel in Sections II A 



Compared to the general definition (see e.g. [3]), two important restrictions appear in (1): (i) 
the expression is the contribution of a single 4-simplex, i.e. it is the first approximation in the A 
expansion of the group field theory; (ii) the area-angle and angle-angle correlations also needed to 
reconstruct the graviton propagator have been discarded. Both restrictions are consistent with 
the choice of dynamics; in fact, the Barrett-Crane model does not give a consistent semiclassical 
limit on more than one 4-simplex, and furthermore it carries no dynamics for the intertwiners — 
the quantum degrees of freedom associated to angles [3, 5]. Therefore, sensible results from 
the BC model can be extracted only if one restricts attention to area correlations on a single 
4-simplex. This in turn means only looking at the "diagonal components" huu{x)h w (y) of the 
graviton propagator. For such components, the BC model proves very useful as a non-trivial 
test for a number of features of the whole approach based upon general boundary formulas such 



II. BARRETT-CRANE DYNAMICS 





and II B. 



as (1). 



A. 



A class of models 



In this paper we consider the class of models defined by the following amplitude, 




(2) 
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where {lOj} denotes the lOj-symbol (described further in Section II B) and the integer k param- 
eterizes the choice of triangle weight in the measure. The interest in this class of models comes 
from the fact that it parameterizes some freedom in the choice of the functional integration over 
the B field. This is the analogue of the usual ambiguity in defining the integration measure over 
4-geometries T>g^ u . The spin foam formalism suggests using the measure that gives exactly a 
topologically invariant partition function before imposing the constraints reducing BF theory to 
GR, a choice which amounts to taking k = 2 in (2). Nevertheless, different choices have been 
considered, and it is interesting to study how they affect the physics. For instance, different 
values of k change the divergent properties of the model, as investigated in [30] for the new 
EPR-based models. Here we study how they affect the next-to-leading order behaviour. 

The form (2) is not the most general BC-like model: one can also consider weights associated 
to each tetrahedron, in the form of powers of the dimension of the intertwiner space associated to 
each tetrahedron [31]. These weights affect the convergence of the model [13, 32]. Such factors 
introduce further coupling between the j;'s and complicate both the analytical and numerical 
studies of the kernel, so we do not consider them here. Notice however that, as shown in [29, 33] 
in the 3d case, both ambiguities should not affect the leading behaviour of the correlations W. 



B. The lOj-symbol and its asymptotics 

The lOj-symbol {10j} is a Spin(4) invariant tensor constructed with Clebsch-Gordan coeffi- 
cients. It is the evaluation of a Spin(4) spin network and can be expressed as an integral over 
five copies of SU(2): 

{10j} = (-l)M f l[dg v n^fenMo)- ( 3 ) 

iSU(2)5 ' 

(For more details see, e.g., [4].) In this expression we associate one SU(2) group element g to 
each tetrahedron v of the 4-simplex, and v(l),v'(l) are the two tetrahedra sharing the triangle 
/. Finally 

sin((2j + l)4>) 

Xj{9) = ^— 

smq) 

is the character of the group element g in the j-representation expressed in terms of its class 
angle <fi £ [0, 27r]. Notice the phase factor (— l)^' 2 -" in (3), which has been so far discarded in 
the spin foam graviton literature [2-6, 8]. We will see in Section III A that it simply corresponds 
to a small change in the choice of boundary state. 

We now briefly recall the known results on the asymptotics of the 10j-symbol, referring to 
the literature for details. The only type of large spin behaviour that has been studied in detail 
is the asymptotics as a fixed set of spins is scaled. More precisely, let ji be ten chosen spins and 
consider the value of the lOj-symbol as the spins are adjusted so that the dimensions 2ji + 1 
become X(2ji + 1) for A a positive integer. 1 In [14] it was shown that certain non-degenerate 



As was already shown in the three-dimensional case with the 6j-symbol [36, 37], the structure of the asymptotic 
expansion is simpler when expressed in terms of 2j + 1 rather than j, and it was also observed in [16] that this 
leads to better behaviour. 
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configurations contribute an oscillating Regge term proportional to 

{10j} ~ A- fl / 2 (-l)S«^ cos A(2j, + 1)0, + Kj), (4) 

l 

where k is an integer and the 6i are the exterior dihedral angles of a flat 4-simplex with areas 
A\ = 2ji + l. We have suppressed a sum over such 4-simplices if there is more than one possibility. 

The geometric saddle point (4) of the lOj'-symbol is notoriously masked by a non-oscillating 
term, which comes from degenerate saddle points in (3) with no direct geometrical meaning. This 
term was found to scale like 1/A 2 and thus dominates the large spin limit of the lOj'-symbol [15— 
17]. The situation is illustrated in Figure 1, where the values of the regular lOj-symbol (with 
2ji + 1 = A) are plotted as a function of A. 



lOj asymptotics, all spins equal 




FIG. 1: Log-log plot of the values of the equilateral lOj-symbol. The dots are the numerical evaluations, 
the lines the asymptotic fits. The oscillating behaviour depending upon the Regge action (predicted to 
scale as 1/A 9 / 2 ) is completely masked by the degenerate configurations, whose leading order is approx- 
imately 2.73/A 2 . After subtraction of this leading order, further non-oscillating contributions appear, 
which are approximately — 5.20/A 3 and 33. 1/A . 

It is interesting to note that even after removing the dominating contribution of order 1/A 2 , 
there are additional non-oscillating contributions still hiding the Regge term. To see this, in 
Fig. 1 we plot the next-to- leading and next-to-next-to-leading orders, obtained by subtracting 
a numerical fit of the previous order. These additional contributions may correspond to those 
predicted in [15], or they may be the power series corrections of the degenerate saddle point 
giving the leading order in 1/A 2 . At the present time, our calculations are not accurate enough to 
be sure that there aren't additional contributions involving fractional powers of A (as predicted 
by [15]), or to be able to discern the expected oscillating contribution of order 1/A 9 / 2 . We hope 
to address this in the near future. 

This situation in which degenerate configurations hide the Regge term prevents a direct 
numerical confirmation of the asymptotics (4), and it has often been seen as a problem with the 
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Barrett-Crane model. The development of the general boundary formalism and the definition 
of the graviton propagator as in (1) has finally allowed us to circumvent this problem, and the 
numerical results we present below are an important test of (4) and in particular of its non- 
trivial dependence on the Regge action. These results are relevant for all spin foam models 
having degenerate configurations emerging in the semiclassical limit, even if those contributions 
hide the geometric configurations, as is also the case for the new EPR-FK-LS spin foam models. 
Following the example of the Barrett-Crane model presented here, the degenerate configurations 
might not be relevant in the study of semiclassical physics with the general boundary formalism. 



III. NUMERICAL STUDY OF AREA CORRELATIONS 

An advantage of the BC vertex amplitude is that its integral representation (3) can be written 
as an integral over ten angles [15, 17], 

{ioj} = (-i)M f n#/c(<« n (5) 
j i i 

Notice that this expression is almost factorized: the ten angles couple only through the (non- 
linear) measure term C(<j)[) = 5(det G[4>[]), where G[4>i] is the Gram matrix of the 4-simplex 
(see, e.g., [4] for more details). This will play an important role in the following. 



A. On the boundary state 

At this point we decide to work with the simplest ansatz for the boundary state, namely a 
factorized boundary state 

^M = T[m, (6) 

where would be peaked around some background value q = (Ao(jo),6), where 6 is the 
dihedral angle. One could also take in (6) different states for each representation label ji, but 
here we make the choice of a homogeneous ansatz with all the ipi's equal, which corresponds to 
the quantization of an equilateral 4-simplex. The factorization allows a simpler analytical and 
numerical study of the correlation. Indeed, using this particular choice of state and the integral 
expression (5) of the kernel, we can rewrite the area correlations (1) as 

w ab (k,a,j ) = ^ j n^cw^^a)^) n ™w> ( ? ) 

I ly^a,b 

where the normalization J\f is j Y\ l dcfti C[<pi\ Y\ l w(<j)[) and the integrands are 

W (4>) = x>j + i) k mxjM, W&) = + ^HtimxM- (8) 

We see that the factorization assumption (6) allows us to perform each sum over the spins 
separately. This turns out to be crucial for improving numerical stability and allowing us to get 
accurate results. 

Concerning the boundary state (6), we consider in this paper the two choices that have so 
far appeared in the literature: 
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• the Gaussian state, introduced in [2] and adapted to the factorized form in [7], where each 
factor in (6) is given by 

• the Bessel-based state introduced in [4, 29], where each factor in (6) is given by 

MJ) = V° l[a) WlUj cos((2j + 1)9). (10) 
V // o(^) " W(^) 

Here \j — jo\ S N and I n {z) is a modified Bessel function of the first kind. 

As we review below in Section VI, the Bessel part of (10) reduces in the large spin limit precisely 
to the Gaussian in (9). Therefore, at leading order, the only difference between (9) and (10) is 
in the phase, which is complex in the first case and real in the second case. Let us recall that 
the rationale for choosing in [4] the form (10) with its real phase was to be able to perform 
the sums (8) analytically, thus turning the definition of the area correlator into an exact group 
integral. This allowed us to obtain an integral representation of the correlator and define the 
perturbative expansion as a saddle point approximation. With the original choice (9), one has 
to first approximate the sums as integrals, and then perform a saddle point analysis of the 
integrals. In both cases, it is ultimately the existence of saddle points with a geometric meaning 
that leads to the desired leading order result. 2 

Here we are interested in a numerical evaluation and thus there is no preference for (10), 
as the sums (8) can be performed numerically for both choices. We first focus on the simpler 
Gaussian boundary state. We discuss the Bessel state and compare the two choices later in 
Section VI. 3 

However, it is important to stress that the assumption of a factorized boundary state is most 
likely too simple to capture a true physical state. The example of linearized Regge calculus shows 
that in general the state will have the same structure of the kernel [40], as already suggested 
in the original paper [2]. However, this factorized boundary state is the only setting in which 
we are able to perform numerical simulations (for large length scales). It is thus a technical 
limitation, and we stress that this is one of the weak points of the present numerical power, 
which deserves attention for the future. 

B. On the numerical methods 

The reason the factorization allows the numerical simulations to be performed much faster 
is that one can bring the sum over the ten spins in (1) inside the integral for the kernel K, 



2 For these geometric saddle points to exist, it is important that the interior dihedral angle 8 — arccos(l/4) of 
the regular 4-simplex be used in the Gaussian state (9). This is due to the extra phase in (3) and (4). Some 
past work [2, 3, 5, 6, 8] used the exterior angle arccos( — 1/4), but those sources also omitted the phase in 
front of (3), and one can show that these two effects cancel out at the leading order of the correlator (1) (up 
to complex conjugation). On the other hand, the phase (— 1)^' 2ji is irrelevant if the spins all have the same 
parity. This is the situation if one uses the Bessel-based boundary state (10), see [4]. In this case, both choices 
of angles in the boundary state give the right asymptotics. More on this below in Section VI. 

3 One advantage of the Bessel choice is to put the integral formulation on more solid grounds, in particular 
allowing one to prove that the degenerate configurations do not contribute to the perturbative expansion [4]. 
A second one is to make the spin ji and the dihedral angle 9i into conjugate variables in the sense of SU(2) 
harmonic analysis. 
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and then decouple them into a product of independent sums. This decoupling simplifies the 
numerical task by decreasing the number of sums over j required to compute the correlations. 
The resulting 9-dimensional integrals are still very difficult to compute numerically, due to their 
oscillatory nature, but can be estimated using Monte Carlo methods with between 10 9 and 10 10 
sample points. This method is able to handle moderately high jo, with sufficient accuracy to 
determine the leading order behaviour. This was used in [7] to support the analytic calculations 
that appeared in the literature, and to give a first glance at the short scale picture. 

Here we push the analysis further, and discuss the next-to-leading order and its dependence 
on the triangle amplitudes in (2). Unfortunately, the Monte Carlo method converges too slowly 
to efficiently extract the next-to-leading order behaviour. An alternative method must be used 
to probe these subdominant asymptotics, which are presented in the next section. This alter- 
native method uses the spin summation formula for the 10j-symbol [34] instead of the integral 
expressions (3) and (5). The result is an expression for the correlations W with no integrals, 
but a large number of sums. Until recently, this spin summation method has been intractable 
due to the large number of operations needed to evaluate the sums. However, the techniques 
developed in [35] have made it feasible and even advantageous to the Monte Carlo method in 
a certain range of jo. The improvement comes from taking advantage of the factorization of 
the boundary state (6) and of hidden structures in the sum to significantly reduce the number 
of operations needed to evaluate it. A more detailed outline of this efficient spin summation 
method is given in the appendix. 

Both the Monte Carlo and the spin summation methods are approximate, the former requiring 
more samples for better convergence and the latter requiring larger truncations of an otherwise 
infinite sum for the same effect (essentially, restricting the boundary spins to lie within a finite 
interval about jo). The difference is that, for jo for which the spin summation method is feasible, 
its convergence is much sharper than that of the Monte Carlo method. In that range, spin 
summation gives at least 5 decimal places of accuracy compared to 2 or fewer for a Monte Carlo 
run of similar length. Fortunately, the feasibility of the spin summation method extends over a 
large enough range of jo and is of sufficient precision to allow us to probe the next-to-leading 
order behaviour of the correlations W as a function of jo . 

C. Numerical results 

We study the area correlator (7) as a function of the squared distance jo between the two 
areas, with attention to how it is affected by changes in the boundary state, i.e. by changes to 
the parameters a and k. First of all, notice that thanks to the symmetries of the background, 
there are only three independent quantities in (7), which correspond to a and b being (i) the 
same triangle, (ii) two adjacent triangles, (iii) two opposite triangles. Here we report the data 
for case (iii). The plots for the other components are qualitatively the same. 

A brief overview of the numerical studies is presented in Figures 2 and 3, for the cases a = 5 
and a = 0.5, and different choices of k. The plots show that the leading order behaviour 
is proportional to 1/jo, with coefficient depending on a. This is compared with the analytic 
calculations, performed following the technique of [2], finding accurate agreement as already 
reported in [7]. We refer to the latter for details on the analytic calculations. We stress the 
fact that measure terms like the ^-dependent triangle amplitudes do not affect the asymptotic 
leading order, in agreement with [4, 29] and intuitively to be expected, since ambiguities in 
the path integral measure should not affect the leading semi-classical behaviour, but only the 
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FIG. 2: Log-log plots of a component of (7) (opposite triangles) for a = 5 and k = 0, 1, 2, — 1 as a function 
of jo. The interpolated dots are the numerical evaluation, with error bars in some cases too small to see, 
and the lines the analytic calculations of the leading order. 

quantum corrections. Indeed, the triangle amplitudes do change the corrections and affect the 
way the asymptotic regime is reached. For instance, it is reached faster for k = 0, hinting that 
this choice minimizes the quantum corrections. We will look at this in more details in the next 
section. 

The plots also show the short scale picture where the non-perturbative regime is probed. 
As presented in [7], the surprising feature is the existence of a maximum in the correlations, 
followed by an inversion of the scaling behaviour: as we go to smaller jo, we reach a peak and 
then the correlations do not diverge but instead are completely regular for small values of jo- 
The precise position of the peak depends on the value of the (Gaussian) width a and on k. We 
will come back to this point in more detail below in Section V. Finally, a closer look at the 
plots shows also some residual oscillations. These are most likely due to the higher order terms 
in the expansion of the Regge action around the background, which contribute an oscillatory 
behaviour as they appear in complex exponentials (see e.g. [2] for details). 

IV. PERTURB ATIVE EXPANSION 

To investigate more quantitatively how the measure term affects the way the asymptotic 
regime is reached, we fit the leading and next-to-leading orders of the numerical evaluations. This 
is a place where the techniques from [35], adapted to the present situation, have been crucial. 
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10° 10 1 10 2 10° 10 1 io 2 



FIG. 3: Log-log plots of the same component for a = 0.5 and k = 0, 1, 2, —1 as a function of jo- 
Computing the graviton propagator using summations over spins gives sufficient precision to 
examine next-to-leading order asymptotic behaviour, unlike the purely numerical (Monte-Carlo) 
integration over group variables. 

After fitting power laws and numerical coefficients, we find that the correlations have the 
following structure: 

W ab {k, a, jo) = - 2S —. V-^-. 2 ho(j ). (11) 

Jo Jo 

Some explicit fits are reported in Table 1. 



(a,k) 


(5,-1) 


(5, 0) 


(5, 1) 


(5, 2) 


(0.5, 0) 


1 /opp 1 


0.127 


0.130 


0.126 


0.120 


0.982 


\ f W | 


0.249 


0.149 


0.241 


0.499 


0.960 



TABLE I: Numerical fits of the LO and NLO for some values of a and k. Only the cases when a and b 
are opposite triangles are reported. 

After the leading order (LO) in 1/jo, the next-to-leading order (NLO) is in 1/jg- This 
linear power increase is natural given the asymptotic nature of the perturbative expansion. 
The NLO in large jo can in fact also be computed analytically, using for instance a saddle 
point approximation as in [4]. However, calculating higher order corrections becomes quickly 



11 



very lengthy and cumbersome, due to the high number of different contributions, already for a 
single 4-simplex. To get an idea, compare the complexity already reached by analogous three- 
dimensional calculations on a single tetrahedron [36], and even just for the underlying 6j-symbol 
[37]. Numerical methods allow us to focus directly on the global structure of the perturbative 
expansion. 

The fitted values for /W reported in the Table approximate well the analytically computed 
values [7], 

/(!)(« = 5) = 0.133, / (1) (a = 0.5) = 1.02. (12) 

From the Table we can observe a slight mismatch on the second digit of the for different 
k. This can be reduced by pushing the simulations to higher jo, however this would be com- 
putationally very expensive, and it is already sufficiently reassuring to see the one digit match, 
together with the plots in Figure 2 showing clearly that the predicted line is being approached. 
The slight mismatch between the coefficients in the table also reflects the fact that the speed 
with which we reach the asymptotic regime depends on the value of k. 

Let us now focus attention on the next-to-leading order. The table shows that the NLO 
depends upon a, just like the LO, but more interestingly it shows a significant dependence on 
k. To highlight this dependence, we plot in Figure 4 the values of f^- 2 \a = 5) as a function 
of k. The plot shows a roughly quadratic dependence with a minimum near the value k = 0. 
This result can be compared with the three-dimensional case studied in [29], where a similar 
quadratic dependence of the NLO coefficient was found and shown analytically. 4 It would be 
interesting (and indeed possible, with enough courage) to compute analytically the NLO in order 
to see the exact dependence. 



0.5 




0.4 




0.3 





0.2 








0.1 





k 



-1.0 -0.5 0.0 0.5 1.0 1.5 2.0 



FIG. 4: The blue dots are the values of f^'{a = 5) at different k. The interpolating line is not a fit, it 
is a simple line linking the plotted values. 



4 For the reader famiiiar with the perturbative expansion of the three-dimensionai geometricai correiations in- 
troduced in [33] and further studied in [29, 36], we would like to point out an interesting difference. Recall 
that both in three and four dimensions, the origin of the good behaviour of the leading order comes from the 
emergence of a Regge-like action from the spin foam kernel, the Wigner 6j-symbol in 3d and the lOj-symbol 
in 4d. In [29] , we distinguished two sources of quantum corrections to the leading order: the ones coming from 
higher orders of the Regge path integral, called (i) in the reference, and the ones coming from higher orders of 
the expansion of the 6j'-symbol, called (ii). In particular, in [29] it was found that only (i) contribute to the 
next-to-leading order on a single tetrahedron, whereas (ii) entered at next-to-next-to-leading order. A similar 
situation might arise in the current case, although we are not in a position to make a precise statement since 
there has not been any computation of the higher orders of the expansion of the lOj-symbol. A difference 
that can be pointed out is that the situation could be even richer, due to the degenerate non-Regge-like terms 
discussed previously. 
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What can we read from the form (11) of the perturbative expansion? Recall the standard 
perturbative calculation of the Newton force [38] , 



The precise value of the numerical coefficients a, does not matter for our considerations. The 
term proportional to aq is a classical correction from general relativity, whereas a2 is the first 
quantum correction. The universality (i.e. mass-independence) of the quantum correction allows 
one to define a running Newton's constant 



The result (13) is computed coupling the graviton propagator to static external matter sources. 
At the state of the art of the formalism, we are far from being able to do this: we need to add 
matter, deal with more than a single 4-simplex, plus have a fully dynamical boundary state. 
Nevertheless, it is interesting to look at the different powers of r and attempt a match with (11). 
The leading order behaviour in 1/jg is consistent with the 1/r 2 scaling of the force in (13). Next, 
the NLO in 1/j'q points directly at the first quantum correction in (13), so again we might say 
that we obtain a consistent scaling of the NLO. Of course, the reader can immediately notice 

3/2 

the absence of a term in l/j i n (H)> which means the absence of the classical correction, i.e. 
the term in a\ of (13). The simplest explanation is that one needs to couple matter to see this 
term explicitly, so its absence should not be taken as further evidence of the inappropriateness 
of the low energy dynamics encoded in the Barrett-Crane model; indeed, the same absence could 
occur when considering the new EPR-FK-LS spin foam models if matter is not coupled. 

Since the corrections to the Newton potential are in principle observable, the discussion 
above shows that fixing uniquely the face amplitudes (in the path integral measure) is one of 
the necessary steps to make a spin foam model predictive. 

V. SHORT SCALE BEHAVIOUR AND METRIC RENORMALIZATION 

The numerical simulations not only allow us to confirm the leading order and study the k- 
dependence of the NLO, but also offer us a look at the possible non-perturbative effects of spin 
foam quantum gravity at short scales, albeit in the simplified model here considered. As already 
presented in the short paper [7], the key point is that, instead of keeping the usual inverse power 
behaviour, the correlations at short scale are regular and do not diverge; they reach a peak at 
some multiple of the Planck scale, after which the correlation recovers its standard behaviour in 
the inverse squared distance. The exact position and shape of the peak depend on both a and 
k, see Figs. 2 and 3. The peak moves (very slowly) to larger spins as we increase the value of 
a, cf. Fig. 5 for the cases a = 50 and a = 100: the location of the peak is now at values of jo of 
order 10. 

These plots also show some unexpected oscillations at an intermediate scale. The frequency 
of these oscillations suggests that they are related to higher order terms in the expansion around 
the background of the Regge action in (4), however we are not able at this stage to give a more 
precise explanation. Let us nevertheless point out that such oscillations already appeared in the 
three-dimensional case for intermediate spins (when looking at the NLO) [29]. In that case, it 
concerned residual oscillations which come from the phase of both the spin foam amplitude and 




(13) 




(14) 
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FIG. 5: Log-log plots for a = 50 and a = 100, k — 0, as a function of j . The lines are the analytic 
calculations of the leading order. The (first) peak is at jo ~ o(10). Notice the presence of oscillations at 
intermediate scales. 

the boundary state, and they could be tamed by changing the width of the Gaussian. We do 
not discuss these oscillations further and instead focus on the short scale picture. 

How can we interpret the presence of a peak, and the very unusual behaviour of the cor- 
relations at scale shorter than the peak? In a truly background independent formalism, the 
semiclassical geometry q is a low-energy property of a quantum state |g), which defines the 
boundary state ^[^afe] = (fla&k); where a and b are three-dimensional indices. It is natural to 
expect that the expectation value of the bulk geometry g^ u depends upon the energy scale set 
by the boundary, i.e. jo- The correlations then allow us to reconstruct the effective metric. It 
is indeed natural in a background independent context to reconstruct the space-time geometry 
and metric from the only available information, namely the correlations (see e.g. [42]). Thus we 
distinguish the flat metric r/ and the actual (effective) metric g. Assuming that the correlations 
W a b scale as l/d g (x, y) 2 in terms of the effective metric, the correlation graph gives the evolution 
of g in terms of the scale d v (x, y) ~ a/Jo (remember that jo is the area scale, thus the square of 
the length scale). In the asymptotic regime for large jo, the metric g converges towards the flat 
metric r/ and we have the right semi-classical limit. For small distances d v (x,y), we have strong 
quantum gravity effects and the effective metric g deviates from 77. 

From this point of view, the correlation peak actually provides a dynamical definition of a 
minimal distance: the physical distance d g (x, y), in term of the effective metric, is never smaller 
than jo 11 ™. This minimal scale result can be compared to similar singularity avoidance results 
obtained in the framework of Loop Quantum Cosmology (see e.g. [39]). 

This dynamical minimal scale j™" depends on the choice of boundary state. It indeed 
depends on the precise value of a. Changing the boundary state would modify the bulk geometry 
and the value of this transition scale between semi-classical behaviour and quantum gravity 
regime. However, we hope that the boundary state and the value of the Gaussian width a will 
be fixed by the requirement of working with a physical state [40, 41]. This would lead to a 
precise prediction for the minimal scale 5 j 7im - 

Actually, at scales shorter than the maximal correlation peak at j™™, the physical distance 
d g (x,y) increases as d„(x,y) gets smaller. This is a sign that the metric becomes more and 
more curved as we get to shorter scales. This is consistent with the expected picture that the 
quantum fluctuations of the geometry create high curvature (and possibly virtual black holes) 
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at the near-Planck scale. We recall that a correlation growing as the one we see is a sign of a 
potential which grows with the distance and thus of confining effects. 

In this respect, let us recall that different approaches to quantum gravity have recently seen 
evidence that the structure of space-time becomes effectively two-dimensional at short scales. 6 
In terms of the area correlations considered here, such a result would show up if the 1/j'o large 
scale behaviour of the correlations is replaced at short scales by log jo. This is not precisely what 
we see, but rather a roughly linear growth. This bears some qualitative similarities with possible 
confinement properties of the theory, but of course any such speculation should be postponed 
to more realistic models. In particular, to obtain the exact effective metric, we need the full 
tensorial structure of the propagator, which can only be obtained using the more complete new 
models. 

Nevertheless this is a crucial line of investigations, and our results show that numerically 
studying the geometric correlations is a fundamental tool for assessing the strength of the spin 
foam formalism to address the long-standing questions of quantum gravity. 



To complete the overview of our numerical studies, we report in this section the results 
obtained using the Bessel-based boundary state (10). Although very reasonable in the large 
spin limit, the Gaussian ansatz can not be considered an exact state in the full theory: the 
conjugation between the dihedral angle and the spin is only achieved taking the latter as real 
variables, an assignment that can be justified only for large spins. On the other hand, in the full 
theory the j\ are half-integers and are conjugated to the SU(2) class angles (see for instance (3)) 
in the sense of harmonic analysis. It thus appears natural to look for a boundary state where the 
conjugation between the spins and the dihedral angles is realized through harmonic analysis over 
SU(2). This is achieved using (10). This state was introduced in [29] for the three-dimensional 
theory and extended to the four-dimensional theory in [4] . The large spin approximation of the 
Bessel functions in (10) gives precisely the Gaussian appearing in (9) [29], 



This asymptotics is reached very early and with great accuracy, as can be seen in Figure 6. 
Therefore this choice of boundary state gives the same leading order of the area correlations, 
and also similar corrections and short scale picture, provided we take the same phases. On the 



5 Here, we have been working with the Riemannian theory and a compact gauge group Spin(4). The representa- 
tion labels j are then discrete. Therefore, we already have a notion of minimal length scale given by the smallest 
non-trivial representation j = 1/2. Notice that our notion of dynamical minimal length is different from this. 
In the Lorentzian theory [11], the representation labels become continuous. Assuming that the correlations 
behave similarly in that case, our definition would lead to a non-trivial notion of dynamical minimal length 
scale. Thus, it would very interesting to numerically check the Lorentzian model to see if it has the same short 
scale behaviour for the correlations. 

6 This has been seen in Causal Dynamical Triangulations [43] and in the Exact Renormalization Group approach 
[44] based on the idea of asymptotic safety [45]. See also [46] for possible similar results in the context of 
quantum group models of quantum gravity, and [47] for a tentative extension of these ideas to loop quantum 
gravity. 



VI. THE BESSEL-BASED BOUNDARY STATE 



/r (^)-/ 2jo+1 (^o) 



V JOTT 2j 



(15) 



15 



o.i 



0.6 



0.5 



0.4 



0.3 



0.2 




5 10 15 



20 



25 



FIG. 6: Comparison of the Bessel state (10) (dots) with the Gaussian state (9) (solid lines), for jo = 2 
and jo = 10. 

other hand, there is an advantage if we change the phase to the cosine given in (10), as in this 
case we can perform the sums (8) analytically [4]: restricting \ j — jo\ S N, we have [29] 



Thanks to this expression, one can now evaluate (8) analytically for k = 0. 7 The resulting 
expression for the correlator (7) takes the form 



This was introduced and studied in [4], to which we refer for details. Obtaining an exact 
integral representation for (1) was the initial rationale to introduce the new boundary state in 
[4], precisely with the hope of improving analytic and numerical control of (1). From the analytic 
point of view, the integral representation was instrumental in proving the conjecture of [2] that 
the degenerate terms of the Barrett-Crane asymptotics do not contribute to the leading order. 
The use of this state was also important from the numerical point of view, as numerics for (10) 
where the first ones to be obtained, confirming the analytic calculations of [4], and paving the 
way to the numerical procedure later used for the Gaussian state, and presented in [7] and here. 

Notice the sums over the rji signs, inherited from the sum in (16). It was proved in [4] that 
each rji configuration gives the same LO asymptotics in l/jo- 8 Therefore for the LO the sums 
over rji give an irrelevant redundancy which is cancelled by the denominator. The calculations 
of the LO were confirmed by the numerics, although the plots were not included in [4]. Here 
we fill that gap and present the relevant plots in Figure 7. The left panel shows the numerical 
evaluation with only the coherent configuration with rji = 1 for all /, whereas the right panel 
shows the result after summing over all the possible rji signs. We see that the LO is the same, 



One can also treat in this way all k even cases, exploiting the property (2j + l) 2n x J (0) = ^%X'{<t > )- General- 
izations to k odd are on the other hand less obvious. 

The coefficients are different than when using the Gaussian state with the complex phase: the difference comes 
entirely from the phase, and it is due to the fact that the oscillating part of (10) is the real part of the oscillating 
part of (9). Because the state enters in both numerator and denominator of (7), the result with (10) coincides 
with the ratio of the real parts of the numerator and denominator computed using (9). 




2k sin 2 (4>~V6) 



(16) 



Wab(k,a,jo) = jr ^2 / Yi d< Pi c i^i\ w va{^a)w Vb ((pb) 11 w m {<j>i) 
{m=±} 1 ¥a,b 



(17) 
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as claimed above, however the full non-perturbative evaluation is rather different. In particular, 
the coherent rji configuration is very similar to the Gaussian state calculation, whereas the sum 
over the rji creates interferences and introduces oscillations which make the global picture quite 
different. In particular, the oscillations also swamp the short scale peak. 

a a a a 

10 1 

10" 1 

lo" z 

io" 3 .... 

10 10 10 10 10 10 

FIG. 7: The Besscl plots. Left panel: the coherent signs configuration, namely the single rji = 1, VZ. 
Right panel: the result after summing over the etas. 

To comment on these oscillations, let us look again at (16). From the Gaussian shape of its 
right hand side, it is clear that this Fourier transform peaks the SU(2) class angle eft on both 9 
and 7T — 9, with spread y a /jo . Notice also that because of the restriction on the spins made 
above, now all ji are either integers, or integers shifted by a half, accordingly to jo. Therefore 
taking 9 to be either the internal or external dihedral angle leads to the same asymptotics, 
a very different situation from the complex phase of the Gaussian ansatz (9), where only the 
right choice of angle in the boundary state leads to a good asymptotics in 1/jo- This symmetry 
between the internal and external angle might look appealing at first, but it is in fact the origin 
of the messier non-perturbative behaviour of (17) highlighted by the right panel of Figure 7. 

Let us recall that the kernel (3), being made of sines, is a product of terms each of which 
contains two exponentials with opposite signs. This is a characteristic of the spin foam formalism, 
which appears also in the new models, and it is interpreted as the fact that the kernel sums 
over orientations of the spacetime manifold. When acting with a state like (9), only one sign 
for each term contributes to the asymptotics, whereas the terms with the "wrong" signs are 
exponentially suppressed. This mechanism shows that a complex phase in the boundary is key 
to give an (observer-induced) orientation to the kernel. Our results show also that fixing this 
orientation is important for having a nice short scale picture without large oscillations. 

VII. CONCLUSION AND OUTLOOK 

In this paper we numerically studied various aspects of the spin foam graviton propagator 
on a single 4-simplex, using the Barrett-Crane model for the dynamics. Although now replaced 
by the better-behaved EPR-FK-LS spin foam model, the Barrett-Crane model still provides an 
excellent model for testing many ideas, and furthermore it is already at the edge of present 
computational power. 

Our results show that numerical techniques are very useful for studying the type of correlation 
described by the propagator. One can go beyond merely checking the analytic calculation of 
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the leading order, and study higher order corrections and short scale physics. We showed how 
both are affected by the triangle amplitudes, corresponding to measure choices in the path 
integral. In particular, the trivial triangle amplitude minimizes the quantum corrections. The 
short scale picture presents the remarkable feature of a peak, followed by an inversed behaviour 
of the correlation. This result, already presented in [7], is very suggestive that interesting 
physical effects can be studied by looking at these types of quantities. We showed how the peak 
depends on the parameter in the boundary state, and then compared the behaviour of this model 
with more realistic short scale scenarios arising from alternative approaches to non-perturbative 
quantum gravity. Finally, we compared the usual Gaussian ansatz for the boundary state with 
the Bessel-based state. 

As a final remark, our results also indicate that it is important to pour effort into overcoming 
the current limitations of numerical simulations, in particular to be able to treat non-factorized 
boundary states, as well as many 4-simplices. 

APPENDIX A: SPIN SUMMATION METHOD 

The spin summation method introduced in Section III and used in Section IV is a general- 
ization of the Christensen-Egan (CE) algorithm for efficient evaluation of the lOj-symbol [34]. 
This generalization is described in full detail 9 , in [35]. Section 3.3 of [35] is most relevant to this 
paper, as it treats the case of a single Barrett-Crane vertex contracted with a factored bound- 
ary state. Below we give a brief overview of these techniques, emphasizing their innovations 
compared to previously known ones. 

To start, we summarize the CE algorithm. As discussed in the original paper, the lOj-symbol 
can be written as a contraction of SU(2) Clebsch-Gordan (CG) coefficients. Each CG coefficient 
is a finite-dimensional rank 3 SU(2) tensor. As such, one can naively try to evaluate the lOj- 
symbol by multiplying all the CG tensor components and then summing over each contracted 
index. Unfortunately, the sums obtained in this way are prohibitively large. The CE algorithm 
makes two significant advances to make the evaluation of these sums feasible. 

First, it uses SU(2) recoupling identities to eliminate the SU(2)-representation index sums in 
exchange for a smaller number of intermediate irrep label sums (spin sums). At this stage, the 
CG tensors have been eliminated in favor of SU(2) 0- and 6j-symbols (see [34] for details), func- 
tions of 3 and 6 spins respectively. These, in turn, can be formally seen as infinite-dimensional 10 
tensors of rank 3 and 6, even though their indices (spins) do not necessarily label basis vectors 
in any linear space. To reiterate briefly, SU(2) recoupling identities are used to express the 
contraction of many CG tensors, with a large number of dummy indices, as a contraction of 
some different tensors, now with a smaller number of dummy indices. 

Notice that some tensor contractions can be evaluated more efficiently than others. In par- 
ticular, the trace of a product of k n x n matrices is a contraction of k n-dimensional, rank 2 
tensors, involving k dummy indices. A naive estimate for the number of operations needed to 
evaluate this contraction as a sum over each dummy index gives 0(n k ). On the other hand, 
standard matrix multiplication algorithms can evaluate the same product trace in 0(/cn 3 ) steps, 



There is an important difference in notation between [35] and this paper. The former uses twice-spins (taking 
on only integer values) to label irreducible representations of SU(2), while the latter uses spins (taking on also 
half-integral values). 

Even though these tensors are infinite-dimensional, sufficiently many of their components vanish to make the 
relevant sums finite. 
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a significant gain in efficiency for k > 3. The second significant advance of the CE algorithm is to 
recognize that part of the expression for the lOj-symbol as a contraction of 9- and 6j-symbols has 
precisely the form of the trace of a product of matrices, which can then be evaluated efficiently. 

Inserting the CE spin sum expression for the Barrett-Crane lOj-symbol into (1) gives a spin 
sum expression for the area correlations W. Thinking of summed spins as contracted tensor 
indices, this expression can be thought of as another large tensor contraction. At first glance, 
its evaluation would require a prohibitively large number of operations, however it is amenable 
to the same kind of analysis as in the CE algorithm. 

Again, two steps are essential in making the spin summation method feasible. First, one 
notices that this new tensor contraction still has the matrix product trace structure of the 
original CE algorithm, however only if the boundary state is factorable, as in (6). This condition 
is where the technical restriction to factorable boundary states appears. Second, the normally 
dense matrices in this product trace can be decomposed into sparse, structured factors. This 
factorization is achieved using SU(2) recoupling identities involving products of 6j-symbols. The 
reason this factorization speeds up the calculations is that multiplication of sparse, structured 
matrices (diagonal, banded and other types) is significantly faster than that of dense ones. 
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